home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / fractal / kaos.lha / complib / linsol.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-11-18  |  2.0 KB  |  98 lines
  1. /*
  2. ### Solving a linear equation, A x = b, ###
  3. -----------------------------------------------------------------------
  4.  
  5.             LINSOL.C Kim 9/1988
  6.  
  7.  
  8.         Algorithm: Guassian elimination with partial pivoting
  9.  
  10.  
  11.         Input: A = n x n matrix, b = n dim vector
  12.         Output: x = n dim solution 
  13.         Return: -1 in case of error, 0 otherwise     
  14.  
  15. -----------------------------------------------------------------------
  16. Bug: Need to specify Maximum matrix dimension.
  17. -----------------------------------------------------------------------
  18. Reference: P136 John R. Rice
  19. ----------------------------------------------------------------------
  20. */
  21. #include <math.h>
  22. int
  23. linsol(x,a,b,n)
  24. int n;
  25. double *x,**a,*b;
  26. {
  27.     int stop,i,j,k,imax;
  28.     double det;
  29.     double r,m,sum,amax,atmp,btmp,fabs();
  30.     static double buffer = { 1.5e-16 };
  31.     
  32.     switch(n){
  33.         case 1:
  34.             if(a[0][0]==0)
  35.                 return(-1);
  36.             else
  37.                 x[0] = 1./ a[0][0];
  38.             break;
  39.         case 2:
  40.             det = a[0][0] * a[1][1] - a[0][1] * a[1][0];
  41.             if(det == 0)
  42.                 return(-1);
  43.             x[0] = (a[1][1] * b[0] - a[0][1] * b[1])/det;
  44.             x[1] = (a[0][0] * b[1] - a[1][0] * b[0])/det;
  45.             break;
  46.         default:
  47.             for(k=0;k<=n-2;k++){
  48.                 imax=k;
  49.                 amax = fabs(a[k][k]);
  50.                 for(i=k;i<=n-1;i++){
  51.                     if(fabs(a[i][k])>= amax) {
  52.                         imax=i;
  53.                         amax=a[i][k];
  54.                     }
  55.                 }
  56.                 if(fabs(a[imax][k]) <= buffer)
  57.                     break;
  58.                 for(j=0;j<=n-1;j++){
  59.                     atmp = a[k][j];
  60.                     a[k][j] = a[imax][j];
  61.                     a[imax][j]=atmp;
  62.                 }
  63.                 btmp = b[k];
  64.                 b[k] = b[imax];
  65.                 b[imax] = btmp;
  66.                 for(i=k+1;i<=n-1;i++){
  67.                     m = a[i][k] /= a[k][k];
  68.                     for(j=k+1;j<=n-1;j++){
  69.                         a[i][j] -= m * a[k][j];
  70.                     }
  71.                     b[i] -= m * b[k];
  72.                 }
  73.             }
  74.             for(i=n-1;i>=0;i--) {
  75.                 sum=0;
  76.                 for(j=i+1;j<=n-1;j++)
  77.                     sum += a[i][j] * x[j];
  78.                 r = b[i] - sum;
  79.                 if(fabs(a[i][i]) > buffer) {
  80.                     x[i] = r/a[i][i];
  81.                 }
  82.                 else if(fabs(r) <= buffer) {
  83.                     x[i] = 0;
  84.                     printf("A is singular, consistent system\n");
  85.                     stop = 1;
  86.                     return(stop);
  87.                 }
  88.                 else {
  89.                     printf("A is singular, inconsistent system\n");
  90.                     stop = 1;
  91.                     return(stop);
  92.                 }    
  93.             }
  94.             break;    
  95.     }
  96.     return(stop);
  97. }
  98.